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ABSTRACT 

The computer program tor the analysis of linearly ela.tlc plane- 
rs. or plane-atrain problems devised by Fellppa In his work oh 
"Refined Finite Element Analysis of Linear and Nonlinear Two-dimensional 
Structures" has been modified to include the use of initial displacement 
boundary conditions. In addition the original IBM 7094 computer dependent 
program has been adapted for use on the IBM 360/65 computer. In both 

programs the FORTRAN IV language has been used. 

Problems involving "Poor fit" displacement boundary conditions and 
refined mesh analysis using coarse mesh analysis input displacements, 
which could not have been done with the original program, are now pos- 
sible with the modified version presented herein. 
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Chapter I 



Introduction 



In this chapter reasons are given for the use of computer structural 
analysis. The purpose of the presentation is explained, and some expected 
program uses are listed. The objectives and extent of work completed by 
the author are presented. 

1. 1 Reasons for computer analysis 

The engineer of the present day is faced with solving structural 
problems of great complexity. Classical mathematics, despite its ever 
increasing sophistication, is only capable of solving severely idealized 
situations while at the same time placing a burden on skilled manpower 
which could be better used for design and development processes. 

Fortunately, the simultaneous development of the digital computer and 
new general methods of solution for problems in Continuum Mechanics has 
come to the rescue in many areas of investigation. The modern digital com 
puter coupled with the finite element method is revolutionizing the ap- 
proach to the process of analysis. Expensive and time consuming experi- 
mental models, now often used in the design of important structures, are 
rapidly becoming displaced by more economical computation. 

1. 2 Purpose of presentation 

This presentation offers a computer program which takes advantage of 
the modern computer's high computational speed and the versatility of the 
finite element method. The program is for the general analysis of arbi*- 
trary plane-stress or plane-strain structural mechanics problems in linear 
elasticity. Sufficient background material on both the finite element 
method and the program structure is given so that those not conversant 
with these matters may use the program effectively. In addition, informa- 
tion on the practical aspects of the program alogrithm for the direct stif 
fness procedures, boundary condition application, and solution of large 
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systems of linear equations are given for those interested in making 
modifications or additions to the program, 

1.3 Expected uses 

Two-dimentional elastic analysis by the finite element method is 
well suited for the development of large scale, production usage computer 
programs. Expected areas of use for the program presented include: 

(1) Practical working design situations, ranging from one-time 
preliminary analyses to extensive in-depth studies. 

(2) Design data generation for stress concentration factors or 
design tables for specific structural shapes. 

(3) Augmentation of other analysis methods. Localized area 
analysis. Utilization of the thermal stress and body force 
features. Gross analysis to determine instrument or sensor 
placement. 

1.4 Scope and objectives 

The original version of the program was written by Dr. Carlos A. 
Felippa [3]* during the course of his Ph.D. studies at the University of 
California. The original program was dependent on the IBM 7094 computer 
and associated hardware then at the Berkeley Computer Center. 

The extent of the author's objectives and work on the program has 
been: (1) reprogramming modifications necessary for the adaptation of the 
program to the faster IBM 360/65 computer at the Naval Postgraduate School 
Computer Center, (2) substituting the use of the more expedient random ac- 
cess disk unit for the original tape storage methods which were used for 

* 

Numbers in brackets refer to the list of references on page 5^ 
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portions of the external storage requirements, (3) restructuring and 
modifying the program to allow usage of initial displacement boundary 
conditions, thus allowing the program to accept all physical boundary 
conditions for planar problems, (4) creation of a "User's Manual 11 to 
facilitate use of the program by persons not familiar with the finite 
element method or production programming techniques* This presentation 
is orientated toward this final objective. 



13 



Chapter II 



General Program Features 



In this brief chapter the general program capabilities and features 
are given as background information. 



2.1 Capabilities 

The program in its current form represents a large scale, high speed 
general computational processor for the stress analysis of plane elastic 
bodies. It is one of the few such programs available in the open litera- 
ture . 



The program is capable of analyzing any arbitrary plane structural 
shape, singly or multiply connected. The body thickness may vary, to a 
degree. All physical boundary conditions are acceptable. Structural 
loading types may include surface forces (concentrated and distributed), 
body forces, and thermal loading effects. Any linearly elastic, homogen- 
ous, isotropic material may be analysed. Up to six different materials 
are allowable for composite material structures. 

2.2 Analytic results 

The method of analysis is the finite element procedure [2], [16], 
with the displacement solution method and direct stiffness matrix genera- 
tion [12], [13]. Displacement compatible finite elements are used. This 
type of element has the virtue in theory that if finer and finer mesh sub- 
divisions are used and there is no numerical error or round off, conver- 
gence to the exact solution is assured [11]. 

In the analysis the program computes in-plane deflections and stresses 

at selected sites on the body, resulting from in-plane loading. Stress 

components (T % O' * and are computed, as well as the principal 

x y xy 

stresses and directions at every nodal point. Stress contour graphs of 
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the structure, or 
structure outline 



a magnified subportion, are generated, wherein the 
and constant stress level contour lines are printed. 
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Chapter III 



Fundamentals of the Finite Element Method 



In this chapter, the fundamentals of the finite element method are 
given in order to acquaint the user with the procedures and the nomen- 
clature needed to understand the program structure and fully utilize its 
potentials . 



3. 1 Basic idealization 

The basic concept of the plane finite element method is that any 
continuous two-dimensional body may be separated by imaginary lines into 
a finite number of individual elements, i.e., an assemblage of smaller 
individual plates. Thus an actual continuous structure is replaced by 
an assemblage of discrete structural elements (Fig. 3.1). 




Figure 3.1 The 
The elements are assumed to be 



Node 




Idealized structure and 
equivalent loads 

Finite Element Idealization 
interconnected only at a discrete 



number of joints or nodal points through which forces are transmitted, 
but the structural elements are such that the connectivity of the struc- 



tural system is preserved along all the common boundries of adjacent 
elements . 

With these stipulations: (1) continuum approximation by an assem- 
blage of elements and (2) element displacement field restrictions, we 
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have, in effect, converted a continuous body into a body composed of many 
structural elements where each element can deform in only a certain number 
of predictable shapes. 

A given point in the continuum has two independent degrees of freedom 
of displacement. The corresponding point in the idealized structure still 
has two degrees of freedom, but the displacement field is now restricted 
to the one selected and imposed on the element. This reduction is basic 
to any solution. In the finite element method the reduction is made in a 
physical manner using structural elements that can retain many of the 
physical properties and behavior characteristics of the original material. 
In other methods of analysis the reduction is often made in a mathematical 
procedure, in which many of the desirable properties and characteristics 
are lost through inability to make adequate or tractable mathematical 
models . 

3 . 2 Displacement fields 

The constraints we have placed on the element take the form of dis- 
placement functions, or fields, for which the overall element displace- 
ment shapes, and resulting strains and stresses, are functionally re- 
lated to the nodal point displacements (generalized coordinates). The 
success of the element in duplicating the performance of the actual con- 
tinuous body is critically dependent on the choice of the displacement 
function. Much effort has been given to the determination of successful 
displacement functions and the fundamental rules governing their genera- 
tion. For flat elements, general procedures have been developed [3], how- 
ever as yet there is no such general method for the construction of dis- 
placement fields for a curved element [1]. 
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3 . 3 Convergence requirements 



To guarantee that a finite element solution will converge to the 
true solution, under conditions of increasingly finer element division 
(finer mesh), the planar displacement function must be able to satisfy 
three basic conditions: 

(1) Continuity condition: interelement displacement compatability 
must be satisfied. For plane cases this implies that adjacent 
element edges must remain together under element deformation. 

(2) Completeness conditions: (a) rigid body modes - overall 
rigid body displacement of the element must not result in 
element straining; (b) constant strain states - nodal displace- 
ments that are indicative of constant strain must result in con- 
stant strain conditions in an element. 

(3) Invariance - element stiffness properties derived through 
the use of displacement functions must remain the same for all 
coordinate systems which are used for their derivation. 

3.4 Element equilibrium equations and stiffnesses 

Knowing element displacement modes as functions of nodal displacements 
we can now write the equilibrium equations for the element. The relations 
between nodal displacements and nodal forces (generalized forces) are ex- 
pressed in compact form by the element stiffness matrix. The resulting 
element equilibrium equations take the form: 



Where , 
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e 




(3.1) 




= Nodal point force vector for an element, (2N x 1) 

= Nodal point displacement vector for an element, (2N x 1) 
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[k] = Element stiffness matrix, symmetric, (2N x 2N) . 

(N is the number of nodes per element) 

The element stiffness matrix is a function of the geometric and physical 

properties of the element. 

The derivation of element stiffnesses can be approached using the 
principle of virtual displacement [13], [16] . Here, virtual nodal dis- 
placements in a local frame of reference are imposed on the element. The 
external and internal work done by the various forces and stresses during 
displacement are equated. The resulting equation can be reduced to the 




[B] = Strain-Displacement relationships from 
elasticity considerations. 

[D] = Stress-Strain relationships from elementary 
theory. 

The element stiffness is developed by performing the indicated matrix 
integration. 

An alternate method of deriving the element stiffness is an energy 
approach [3], [5]. Here the total potential energy of the element-load 
system is minimized, an outgrowth of the application of the governing 
variational principle. This is, in fact, a form of the Ritz technique 
applied to the network of finite elements. This method is not as direct, 
nor as physically interpre table as the previous approach. However it is 
a more powerful technique and leads naturally to the formulation of the 
consistent mass matrix, consistent stiffness matrix and consistent load 
vectors. Elements having sophisticated geometries and higher order dis- 
placement fields can be handled practically only with the energy method. 
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Having the element stiffness matrix and the two auxiliary matrices 
[B] and [ D ] we can develop the full stress analysis for a single element 
when the nodal displacements resulting from the loading are known. 

3.5 Complete stiffness matrix 

To permit combining the individual stiffnesses each must first be 
transformed from its local coordinate system to a common (global) system 
for the entire node-element mesh. This is done with a conventional tensor 
transformation relating the two coordinate systems. The formulation of 
the complete stiffness matrix for the discretized structure can now be 
accomplished. In the direct stiffness method this is achieved by the 
direct addition of the element stiffnesses at all the element interface 
nodal points. This is a technique idealy suited for computer operation 
(see appendix 2). 

3 . 6 System equilibrium equations 

The complete stiffness matrix [K] is the matrix of equilibrium equa- 
tions for the total idealized structure. The set of overall equilibrium 
equations is of the same form as equation (3.1), viz.: 




(3.3) 



The nodal point force vector (load vector) contains the loads to which the 
body is subjected. The loads are in the form of concentrated nodal point 
forces directed along the global coordinates axes. Where the body loading 
is in the form of distributed forces, body force loads or thermal loads, 
equivalent concentrated forces are applied at the appropriate nodal points. 

The nodal point displacement vector (displacement vector) is the un- 
known in equation (3.3). This is a direct result of the formulation 
developed and is responsible for the name, ‘displacement Solution Method.'* 
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Another general method of finite analysis is available, where a stress 
field is used [5], this procedure is called the "Equilibrium Method. 11 

The overall equilibrium equations must be modified by the boundary 
conditions of the problem* This involves elimination of stiffness contri- 
butions in [K] for constrained nodes. In the case of initial displace- 
ment boundary conditions the load vector and displacement vector must be 
adjusted to reflect equivalent nodal forces and the initial displacement 
value, respectively, that result from displacing a node a known amount. 

The final set of equations (3.3) can attain rather large proportions. 
The order of the load and displacement vectors and the size of the sym- 
metric overall stiffness matrix are 2 x P, where P is the total number of 
nodal points in the mesh. Fortunately the stiffness matrix is sparsely 
populated and banded, which permits specialized solution techniques on a 
computer. Actually, the displacement solution method and the direct stif- 
fness procedures are used so that these stiffness matrix properties may 
be exploited (see appendix 2). 

The final equations (3.3) are solved for the nodal point displace- 
ments. With the known displacements, the matrices [ B ] and [D] may be 
applied for each element and the overall nodal stresses evaluated. The 
stress levels at nodes common to more than one element are usually averaged 
to obtain a representative value. The end result is the displacements 
and stress levels at each node of the idealized structure. 
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Chapter IV Program Information and Structure 



This chapter offers information on the program and gives its 
structure, functional routines and general method of operation. 

4. 1 Program identification 

PSELST - Plane Stress Elastic Analysis using Linear Strain Triangles 
Programmed: Carlos A. Felippa, June 1966 

Revised: John P. Malone, July 1968 

4.2 Purpose 

The purpose of the program is to provide a high speed, production 
use computational solution of general plane-stress or plane-strain static 

r 

linear elastic problems using linear strain triangles in the finite ele- 
ment method. Surface loads, body forces and thermal effects may be con- 
sidered . 

4.3 Programming information 

The program is written in FORTRAN IV language [6] for the IBM OS/360 
Model 67 computer. An overlay structure [7], [8], under control of the 
system linkage editor, is utilized to conserve main core storage. The 
overlay segments are arranged so that only active subroutines and atten- 
dant internal storage requirements are in main core at a time* 

4.4 External storage 

Fortran logical units 1,2, 3, 7,8, and 9 are used for temporary stor- 
age. Logical unit 7 is defined as a data set on one IBM 2311 random ac- 
cess disk storage unit. The remaining logical units are defined as 
separate data set blocks on a second 2311 disk unit. The input card 
reader, output printer, and output card punch are defined as units 5,6, 
and 4 respectively. 
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4.5 Basic finite element mesh units 



The basic mesh element is a quadrilateral composed of four 6-nodal 
point linear strain triangles (LST) , the center point being the centroid, 
(Fig. 4.1). Internal points 9 to 13 are eliminated by matrix condensation 
[3], [9] thereby reducing the number of degrees of freedom from 26 to 16. 



Single 6-nodal point triangles may also be specified to facilitate 



fitting of certain shapes. 




(4 corner nodes, 8 external nodes) 
(#9 = centroid) 



KEY : 

—O— 



A 

i — ^ 



Corner Nodes 

Midside Nodes 

Removed during 
condensation 




(3 corner nodes, 6 external nodes) 



Figure 4.1 Basic Finite Element Mesh Shapes 
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4.6 Capacity 



The mesh input is subject to the following limitations for an IBM 



computer with 256 K bytes of main core storage; 

Max. number of elements (Quad, and/or triangle) - - 350 
Max. number of external (total) nodal points - - - 1050 
Max. number of restrained components ------- 250 

Max. difference of nodal point numbers for 

the same element ----------------- 79 



These limits are dictated by storage requirements in the computation 
of the overall stiffness matrix, and not by the equation solver. 

It is estimated that with 512 K bytes of main core storage the pro- 
gram could be modified to accept 1700 external nodes and a half-band- 
width of 250. This would result in approximately 600 usable elements. 
Solution times for full scale problems with these modifications may be 
prohibitive on current computers. 

The most commonly encountered limitation of the current program ver- 
sion is the maximum nodal point number difference. Next in order are 
maximum external nodes (total nodes), maximum number of elements, and 
maximum number of restrained components. 

The maximum number of degrees of freedom is 1050 x 2 (external nodes) 
+ 350 x 10(condensed and recovered nodes) = 5600. The condensation pro- 
cedure reduces the number of equations to 2100, yet retains the versatil- 
ity of 5600 degrees of freedom. 

The maximum half-bandwidth (2 x Max. node number difference + 2), a 
measure of the width of the banded area of the overall stiffness matrix, 
is 160. 
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4.7 Program structure 



The overlay root-segment structure is shown in Figure 4.2, where 
each rectangle represents a subroutine. The subroutine functions are: 



MAIN 


remains in core and controls calling sequence during execution 


RDDISK- 

WRDISK 


remains in core and provides high-speed deposit and retrieval 




of data on a random access disk storage unit; 


SETUP 


inputs, prints and checks mesh data, and evaluates element 



stiffnesses ; 



STQUAD 


assembles and condenses quadrilateral stiffness; 


STLST6 


computes stiffness of a six nodal point triangle; 


LDINPT 


inputs load cases and reduces surface, body, and thermal loads 



to equivalent forces on element external nodal points; 



THERLD 


computes initial thermal forces for a single triangle; 


FORMK 


assembles the complete stiffness matrix; 


SOLVE 


obtains nodal point displacements from BIGSOL; 


BIGSOL 


solves large capacity banded matrix problems; 


STRESS 


evaluates and prints element and nodal stresses; 


TRISTR 


computes stresses for a single triangle; 


CNTPLT 


produces printer plots of stress contour lines. 



The subroutine flow chart is presented in Fig. 4.3. 



25 



ttfSIQHM/MSia 



cOi 




LO 




* 



-P & 

0) -P 
a) rd 
a.p. 

o ' ^ 



-P 

O 

O 

0$ 



-P 

c 

<D 



a) 

co 



(N 





CO 

CO 

W 


00 

r- 


Eh 

CO 

M 

Eh 


& 

Eh 

CO 












& 




W 


i-t 


P3 


> 


— ' 


O 






CO 


O 




o 


CO 




M 






CQ 



CN 




<D 

tr 

rd w 

m a) 

O -P 
-p >1 
W ,Q 



I 

O 

iH 

•H 




C C 
•H -H 
rd 

E w 

-P 

rH C 

rd a 

TO O rd 
H p g 
H »H 

3 u 
^ D 1 o 
G <D TO 
p w 



> 

•H 




26 



Figure 4.2 Program Overlay Structure 



- START - 



All transfers 
indicated by- 

are controlled by 
MAIN, and require' 
loading and 7 

unloading of 
overlay segments 




no 






Figure 4.3 Subroutine Flow Chart 
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4 . 8 Programming techniques 



Three programming techniques that are contained in the original 
version of the program are noteworthy: 

(1) RDDISK/WRDISK Subroutine - used as a data transfer device to 
external storage. Benefits are: simplified data address accounting 
in subroutines, fast data transfer, access to large storage areas. 

(2) Matrix Compression (condensation) - upper half-band of the 
overall stiffness matrix is stored, column-wise, by blocks for 
ease of manipulation and to reduce storage requirements. 

(3) Use of one dimensional array - to minimize address calculation 
time. 

Techniques (2) and (3) are described more completely in appendix 2. 
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Chapter V 



User Information 



In this chapter specific features of the program are explored. 
Techniques are given to assist in the effective use of the program. 

$•1 Typical applications 

Representative applications of specific problem types are illustrated 
in Fig. 5.1. Case (a) is a wing frame where gust loading is simulated by 
a static load situation. Thermal loading effects at the frame tip might 
also be considered. Case (b) is a cross-section of a length of pressure 
piping. The circular perimeters have been distorted during a bending 
process. The stress levels and patterns at the thin wall section are 
desired. Case (c) is a comparison study of two welds. The left fillet 
welds are cheaper to produce in some given process, but the right side 
full-penetration welds are expected to be significantly stronger, A 
cost effectiveness study requires information comparing the two welds. 

Case (d) is a stress concentration factor study. Data for design tables 
is desired to guide design engineers confronted with this configuration. 
Loading situations of tension, compression and in-plane bending might be 
considered. Case (e) is a section of a prefabricated tunnel proposed for 
underwater installation. The tunnel is partially evacuated as part of a 
high speed transportation scheme. The optimum tunnel cross section is 
desired, using a minimum material criterion. Case (f) is a hypothetical 
plane-strain problem illustrating the versatility of the finite element 
technique used in conjunction with a high speed computer. 

5.2 General p rogram input and output 

✓ 

The versatility of the program is obtained by accepting a penalty of 
requiring a large amount of input data. The data is in the form of stand- 
ard punched cards. The number of input cards will be between 25 and 500, 
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dependent on the number of elements used and the extent of loading. 

Input information is divided in two groups, structure data and load- 
ing data. The general content of input data follows: 

1. Structure data. 

a. Accounting totals - number of elements, nodes, etc. 

b. Mesh configuration. 

JL. Node-to-element identification. 

2. Corner nodes, X-Y coordinates. 

c. Material physical constants. 

d. Boundary conditions 

2. Loading data. 

a. Accounting totals - number and types of loads. 

b. Load magnitudes and locations. 

Once a problem has been set up with the structure data input, repeated 
load cases may be run. Consecutive problems may also be run. 

Program output is divided into three parts: (1) echo check of input 
data, (2) displacements and stresses, (3) contour maps. The general con- 
tent of each part follows: 

1. Echo check. 

a. Input structure data printed for each problem. 

b. Input loading data printed for each load case. 

2. Displacements and stresses. 

a. Final load vector. 

b. Nodal point displacements printed (punch option). 

c. Computed element nodal stresses (optional). 

d. Averaged nodal stresses, CT , rf > T > CT , O', > 

° x w y xy max rain 

T > and maximum principle stress orientation (punch 
max 

option) . 
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FIGURE 5.1 




Typical Applications 
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3. Contour graphs 

a. A contour graph for each averaged stress computed 
(optional) . 

Details of input, output and program options will be covered in subsequent 
sections. 

5.3 Precautions 

This section is presented here to prevent users from rushing to a 
card punch machine with this manual in their hand. 

The author has found that the following ,, rules n must not be violated 
when using the program. 

(1) Know your problem. 

(2) Be patient. 

(3) Be prepared. 

An example of the first rule involves how much of the structure to 
analyze. Many structures have 1-fold symmetry, less frequently 2-fold. 

If the loading is symmetric about the same axes, partial body analysis 
may be utilized. This avoids duplication of data and allows more ef- 
fective use of the elements available. 

The rule on patience should be applied during mesh construction, 
nodal point numbering, and card punching. The program has many error 
printout exits that determine and illuminate input mistakes, but these 
exits also terminate program execution. 

Being prepared means having a plan to make effective use of the 
tremendous amount of information that the program produces. A cursory 
glance at the sigma-max graph to determine if a material is suitable may 
fail to make adequate use of the wealth of output information available. 
Other persons may need information from your analysis or some unusual 
condition may exist that warrents further study. 
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5.4 Mesh construction 



An effective mesh is one that represents a good approximation to the 
true structure. How effective a mesh is may not be known until one run 
has been completed and stress patterns are reviewed. 

An efficient mesh is one that is formed on the basis of the above re- 
quirements, but also allows computation to proceed as rapidly as possible. 
This section deals with the construction of efficient meshes. 

The mesh idealization of the body may contain 1 to 350 basic mesh 
units (Fig. 4.1). Due to its superior stiffness properties, use of the 
quadrilateral is preferred. The use of triangles should be limited to 
curve fitting difficulties and when changing the fineness (or coarseness) 
of the mesh; this is termed grading. 

The nodes of the basic mesh unit are called corner points or midside 
points. Mesh units may only be joined with correspondance between corner 
and midside points. Two adjoining units share one midside point and two 
corner points. In a mesh, corner points will join three (usually four) or 
more elements when internal to the body (Fig. 5.2c, points a and b). Mid- 
side points are always common to only two elements when within the body 
interior. Each element and node must be numbered for identification pur- 
poses. Element size, shape and location information is determined by 
listing the X-Y (global) coordinates of all corner points. Midside point 
coordinates and element local coordinate systems are internally generated 
in the program. Nodal numbering techniques and coordinate system usage 
are covered in the following sections. 

Three factors should be considered in constructing the mesh: 

(1) Expected or known stress gradients. 

(2) Mesh symmetry. 

(3) Nodal point numbering. 
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Stress gradients should dictate the degree of fineness of the mesh. 

In areas of high stress gradients the mesh should be refined. Portions 
where stress levels are nearly constant can be represented with large 
elements in a coarse mesh. In many problems users intui tative ly produce 
meshes that account for stress gradients. However, when complex shapes, 
boundary conditions and loads are involved, a preliminary coarse mesh 
analysis should be conducted. The indication of high stress gradients is 
not always available from the contour graphs, since they are based on 
averaged values of common nodal stresses. A true indication of the need 
for element refinement is when element stress levels at shared nodes 
differ considerably. 

Many techniques for grading the mesh during initial construction, or 
when refining is required, are available. A comparison of two grading 
methods is shown in Figure 5.4. 

Mesh symmetry with respect to some line or curve in the body should 
be maintained where possible. In complex shapes overall symmetry usually 
cannot be achieved. In these cases local symmetric areas should be con- 
structed. Figure 5.2 illustrates symmetric mesh construction. 

In a completed mesh the one element that has the largest difference 
between its nodal point numbers influences heavily the computation time 
required for the problem. The maximum difference in node numbers is an 
indication of the amount of storage required to house and solve the set 
of equilibrium equations. In solving the equations the data in storage is 
transferred, manipulated, and restored many times. The larger the data 
set, the longer the corresponding computation time. The mesh can be con- 
structed so that it lends itself to optimum nodal point numbering tech- 
nique, which in turn will reduce the largest nodal point number difference 
and the computation time. 
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KEY: 




Unsyminetr ic (poor) Mesh 




Improved Symmetric Mesh 
(reduced computation time) 




Corner Nodal Point 



Midside Nodal Point 




Figure 5.2 Symmetry in Mesh Construction 



35 



In general, the least nodal number difference can be achieved by plac- 
ing all mesh corner nodes on lines or smooth curves within the body. Figure 
5.2c illustrates this procedure. 

It is suggested that the physical process of mesh construction be ac- 
compolished in the following manner. (1) Sketch initial mesh ideas on a 
scratch pad until the final version is formed. (2) Check the limitations 
on maximum numbers of elements, nodes, and constraints and the maximum 
nodal number difference. (3) Transfer the mesh layout to a large sheet 
where sufficient space is available to identify corner and midside nodes 
by number and to exhibit the X-Y coordinates of corner nodal points. Large 
(2 , x3 l ) Ozalid prints from a master drawing of faint \ inch square grid 
with a superimposed quadrant of concentric arcs have been found convenient 
for this last step. 

Meshes can be produced where all the elements are of equal size and 
have the same orientation. In this case all elements have the same stif- 
fness matrix and the program need compute only one. This is a time saving 
feature that can be used for regular bodies or preliminary gross analyses. 
Frequently an auxiliary computer program can be written to generate all 
mesh input data cards in these cases. 

5.5 Nodal point and e lement numbering 

The program requires that each element and each nodal point be numbered 
for identification purposes. The numbering of elements and nodes is in- 
dependent, but increasing both counts in the same pattern is convenient. 
Consecutive numbering is not a requirement, but again, it is convenient. 

If consecutive numbering is not used (occurs when combining two previously 
independent meshes) an ordered listing of the input sequence should be 
maintained, since output is identified by consecutive numbering. 
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The method of nodal point numbering is important for two reasons, 

(1) To keep within the limit on maximum nodal point difference 
for a single element. 

(2) Efficient numbering can significantly reduce computation 
time. 

The general rule for efficient numbering is to number in a direction 
where the least number of nodes lie in '‘line" (not necessarily the small- 
est figure dimension) and repeat for the next ‘‘parallel line." The 
technique is illustrated in figures 5.3 and 5.4c. 




Figure 5.3 Nodal Numbering Technique 
Efficient numbering within graded areas is often difficult. Figure 5.4 
compares two grading techniques and three numbering methods. Figure 5.4(c) 
is preferred on both scores. 

As discussed previously, the one element in the mesh with the greatest 



37 



node number difference controls the solution time. If it is unavoidable 
that a large nodal number difference exists in some area of the mesh, 
nothing is to be gained by trying to reduce lower difference values else- 
where in the mesh. 

5.6 Coordinates 

The global coordinate system for the program is the planar cartesian 
coordinate system. The X and Y coordinates of each comer point are re- 
quired input data. This coordinate array data set is usually the largest 
single block of input data. Two hints are offered that can reduce the 
time to prepare these data cards. (1) Reduce the significant figures of 
the coordinate by slight nodal point movement. This procedure may disrupt 
mesh symmetry, but the effect will be small. (2) The coordinate origin 
can be placed internal to the figure. Mesh symmetry may be utilized to 
produce identical coordinate values, except for sign. 

The coordinate system is compatible with structure overlapping. This 
situation has not been explored other than to demonstrate its dreadful in- 
fluence on the contour graphs, which should not be printed in these cases. 
Overlapping can also be produced as a result of deformation and the pro- 
gram does not sense the condition. 

5.7 Units 

A consistent set of units must be used between the linear measure and 
surface traction load values. Where the inch is used as the unit of bound- 
ary and thickness measure, distributed forces and shears must be per square 
inch (eg. PSI). The modulus of elasticity and specific weight input con- 
stants must reflect the same units. The most useful set of units for 
mechanical design has been the inch and KSI (i.e. Kip per square inch); for 
large structural problems the foot and KSF. 
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Figure 5.4 Grading and Numbering Techniques 



5.8 Materials 



The program is written for homogenous, isotropic, linear elastic 
materials. The required input data are the Elastic Modulus, E; Poisson's 



The ability to utilize multiple materials is one of the advantages 
of the finite element procedure. In the program up to six different 
materials may be specified. True difference between materials may be 
used, or different property values of the same material. Temperature 
influence on the elastic properties may be approximated by using "differ- 
ent" materials when in the thermal gradient of a single material object. 
Body force loading, other than the program standard of 1-G, can be simu- 
lated by applying suitable multipliers to the value of the specific weight 
prior to input. This technique is useful in simulating dynamic inertia 
loading by a static equivalence. 

When multi-material problems are analyzed, the element stress print 
option should be specified. The standard print of averaged nodal point 
stresses is not adequate at material interfaces. The stress contour graph 
should be used only with judgement in these cases, since they do not re- 
flect true stress discontinuities at interfaces. 

5.9 Boundary conditions 

Boundary conditions are applied at the appropriate nodal points 
(corner or midside) of the structure. A boundary constraint consists of 
fixing a node in one direction. A "fixed" boundary point will require two 
boundary constraints. The program will accept up to 250 constraints. 



Ratio, 7/ ; Specific Weight,^ ; and Coefficient of Thermal Expansion, C?C 
For plane-strain problems reduced (primed) values must be used: 
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Constraints are not restricted to nodes on the physical boundary; how- 
ever, this is the usual situation. 

Initial values of boundary displacement may also be used as boundary 
conditions. This feature has many applications. One is the simulation 
of "poor fit 11 conditions in a part or structure, A second is illustrated 
in Figure 5,5, where actual gear tooth root area deflections were deter- 
mined; then used as input initial displacements, thus freeing elements for 
use in a more refined mesh of the structure of interest. 

In program input the boundary condition data is designated for the 
appropriate nodal point with a boundary condition "TAG." A list of TAG 
values and appropriate boundary conditions follow: 



TAG VALUE 



BOUNDARY CONDITION 



TAG = 0 



If the point is fixed in both directions. 



or 



Initial X and Y displacements are specified and 



the point is fixed in both directions. 



TAG = 1 



If the point is fixed in the X-direction and 



free to move in the Y-direction, 



or 



Initial X-displacement is specified but the point 



is free to move in the Y-direction. 



TAG = 2 



If the point is free to move along a line forming 
an angle (p with the positive X-axis. 



or 



Initial Y-displacement is specified but the point 



is free to move in the X-direction. 



The boundary conditions are illustrated in Figure 5.6. 
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Initial coarse mesh 
determines displacements 




Subsequent refined mesh. 




Figure 5.5 Use of Initial Displacement Boundary Conditions. 
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Nodes 1,2,3 TnG=0 
Initial X-DISPL =0.0 
Initial Y-DIS n L =0.0 






Nodes 1,2,3 TAG=1 
Initial X-DISPL = (-d) 



Nodes 1,2,3 TAG=2 
<£>=Bour.dary Angle=150.0° 

5) Initial DISPL = 0.0 (Required) 



Nodes 1,2,3 TAG=2 
<p=Boundary Angle=0 . 0° (Required) 
Initial Y-DISPL = (-d) 









Initial Y-DISPL 






Figure 5.6 Boundary Conditions 
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A single initial displacement can only be specified in the X-direc- 
tion (TAG=1) or the Y-direction (TAG-2, = 0.0°). Initial displacements 

for points on a slope roller are not available in this version of the pro- 
gram. There is an equivalance between TAG^l and TAG=2, (f) = 90.0°. The 
first method requires less internal computation. 

In multi-load case problems any initial displacements specified will 
apply only for the first load case. Subsequent load cases will have the 
same constraint conditions for initially displaced nodes, but with a 0.0 
value of initial displacement. 

Boundary constraints are a fundamental characteristic of a structure 
and are known immediately upon the definition of a problem. When the ad- 
vantages of 2-fold or 4-fold symmetry are utilized, with partial body 
analysis, the constraints needed are not always obvious. Care should be 
exercised in these cases to ensure duplication of the actual problem 
situation. When a complete structure is analyzed that is shape and load 
symmetric, auxiliary boundary constraints are often helpful to improve the 
solution. Frequently symmetric centerline nodes can be placed on rollers. 
This will avoid slight skewing displacements which result from round-off 
error in the stiffness matrix generation* 

5.10 Loading 

Three loading types are available with the program: (1) surface loads, 
(2) body forces, and (3) thermal loading. The general character of the 
input loading data is listed below: 



Loading type 
Concentrated force 

Distributed normal force* 

Distributed shear force* 



Input characteristics 

Magnitude of X and Y components specified. 
Loaded node identified. 

Magnitude and sense of normal pressure. 
Element side identified. 

Magnitude and sense of surface shear. 
Element side identified. 
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Gravity (body force) load Indicated true (T) if applicable. 

Acts in (-Y) direction. 

Thermal load Indicate temperature increment at 

element comer nodes. 

* - side variation assumed parabolic based on values specified 
at the nodes of the element side. 

The sign convention for concentrated forces is that of the global 
coordinate system. Distributed forces follow a sign convention based on 
the element shape. Normal outward traction and counterclockwise shears 
are positive (Fig. 6.2). Gravity loading acts in the (-Y) direction. 

Distributed forces and shears acting on an element side are assumed 
to have a parabolic variation based on the input values at the corner 
nodes and midside node involved. The resulting total element side load 
is represented by equivalent nodal point concentrated forces that are 
generated internally in the program. When the gravity load option is used, 
the program generates nodal forces equivalent to body forces to apply at 
all the mesh nodes. In thermal loading the program computes the forces 
that would be required to reestablish the size of the thermally expanded 
elements and applies them at the appropriate nodes. These computations 
are internal to the program, with no external indication other than the 
resulting displacements and stresses. 

5.11 Contour graphs 

The stress contour plots generated by the program are very useful. 

The plots are produced on a conventional printer, thus eliminating a re- 
quirement for specialized plotting equipment not available at many computer 
installations. The loss in definition due to printer plotting is compen- 
sated for by the short time required to produce the plots (about 9 sec.). 

Six plots are available, each shows the body outline and contour lines 
of constant stress. The lines divide the full range of stresses encountered 
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into ten equal increments. Plots of the following stresses are available; 

^x’ Oy* ^xy* ^max* ^min 5 ^max* 

The plotting subroutine allows up to 50 elements to be eliminated 

from the area considered in determining the stress range for the contour 

lines. This feature allows (1) the plotting of only a portion of the 

structure, (2) elimination of high stress gradient areas that bunch the 

countour lines, (3) a combination of the two. 

5.12 Error exits 

Several error conditions caused by either wrong input data or exceeded 
1 program capacity are checked in subroutine SETUP. Error messages printed 
are self-explanatory and may be complemented by examination of the input 
data printout. The program does not stop until all error conditions have 
been tested. If another problem follows after an error detection, the 
program searches for the next START indicator card at which time execu- 
tion continues on the new problem. 

5.13 Displacement solution iterations 

Double precision iterations for improvement of the displacement solu- 
tion are optionally available in the program. One or two iterations 
should always be used for large ill-conditioned problems. Each iteration 
is essentially re-solving the displacement equations to eliminate round- 
off error, and will take approximately 90% of the original problem solu- 
tion time. 

5.14 Card punching of displacements and stresses 

A program option is the punching of nodal displacements and averaged 
nodal point stresses. These card data sets may be used as input to plot- 
ting equipment for the production of more sophisticated plots than pro- 
duced by the program. 
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5.15 Timing 



Sufficient problems have not been run with the program to develop an 
accurate empirical relationship giving overall execution time. An esti- 
mator currently in use is: 

60 nodes, half -bandwidth < 60 

1 minute execution time for each: 50 nodes, 60 ^ half-bandwidth < 120 

(no iterations) 

40 nodes, 120 ^ half-bandwidth 

For small problems (number of elements < 40) execution times pre- 
dicted will be low since the subroutine SETUP requires a minimum of 50 
sec. regardless of problem size. Timing estimates apply to an IBM OS/360 
computer. 

5.16 Accuracy 

To illustrate the accuracy of this version, a 2:1 symmetric plate 
under uniform in-plane load on top (Fig. 5.7) and constrained by dia- 
phragms at the ends (X-disp. only) was analyzed by subdividing one-half 
of the plate into meshes of (n x n) quadrilaterals - square in this case. 
There is no known analytic solution, but this problem has been used ex- 
tensively in comparison tests of different types of finite elements be- 
cause of the diversity of stress conditions and the ease of preparation. 
Precision is not sought in the analysis, since a ,, true n solution is 
guaranteed if the mesh is sufficiently refined. Rapidity of convergence 
is the desirable characteristic. 

A comparison of the typical values reproduced in Figure 5.7 shows 
that the 8x8 mesh provides almost 5 decimal digits for the displace- 
ments, 4 for CT^ and 3 for T" X y* Actually, the program has capacity for 
an 18 x 18 mesh (5500 degrees of freedom) if necessary. The consistency of 
the stress values is reflected by the fact that the maximum discrepancy over 
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contributing elements 



rf , and 0.006 for 

y 



nodal points did not exceed 

7~ in the case of the 8 
xy 



0.004 for O' and 
x 

x 8 mesh. 
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Y 



p ■ 0.1 ksf 



* 




Mesh of square 
elements for 
one-half of plate 


2x2 


h x k 


8x8 


Asymptotic 

Extrapolated 

Values 


Degrees of Freedom: 
unconstrained 


kZ 


290 


1090 




after B.C.s 


32 


272 


1056 


- 


Deflection x 10^ 


U.8297 


U.8U87 


U.8515 


1u8516 


Normal Stress 0*g 


0.3190 


0.320? 


0.321U 


0.3215 


Shear Stress T , 
xyA 


0.1550 


0.1380 


0.1333 


0.1321 



fig. 5.7 - Example: 2:1 Symmetric Plate 
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Chapter VI 



Input Data Preparation 



This chapter is written as a self contained unit giving all specific 
information for the preparation of input data cards. When facility with 
the use of the program has been attained, this chapter may be duplicated 
and used as a condensed program user's manual. 



FLOATING POINT (F) FORMATS ARE COMPATIBLE WITH ASSIGNED EXPONENT (E) FORMATS 
RIGHT JUSTIFY INTEGER AND EXPONENT NUMBERS IN THEIR ASSIGNED FIELDS 
6 . 1 Structure data 

(a) Start Card (A8) : With the word START punched in cols. 1-5. 

This card must precede the input data deck of any problem. 

(b) Title Card (20A4) : Alphameric information in cols. 1-80 to 

identify the output. 

(c) Control Card (814, 6L2) : 





Variable 




Columns 


Name 


Meaning 


1-4 


NUMEL 


Number of elements 350); 


5-8 


NUMCP 


Number of comer points; 


9-12 


NUMNP 


Number of total nodal points ( >$ 1050); 


13-16 


NUMBC 


Number of restrained nodal points; 


17-20 


NUMPB 


Number of plot defining boundary points 
see item (e) ; 


21-24 


NLOAD 


Number of load cases; will be set=l if 
left blank; 


25-28 


NMAT 


Number of different materials 6); 

will be set=l if left blank; 


29-32 


MAXIT 


Maximum number of residual iterations 
in the displacement solution; punch a 
1 or 2 for large, ill-conditioned prob- 
lems . 


The 


next five fields are 


for logical flags; if a T is punched, in any 


assigned 


column, the indicated 


option takes place. A blank of F implies 



FALSE. 
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T1 All quadrilaterals have the same stiffness 

matrix (see note 1); 

T2 Punching of mesh nodal point coordinates 

and displacements (14, 2F8.3, 2E14.5); 

T3 Punching of averaged and IT ^ 

at mesh nodal points and quadrilateral 
centroids (14, 3E18.6); 

T4 Print of element stresses (see note 2); 

T5 Another complete problem follows* 



(1) All quadrilaterals have the same stiffness if they can be super- 
imposed by a translation. 

(2) Element stresses should be printed in problems involving several 
material types, since averaged stresses and their plots do not 
display actual interface discontinuities. 

(d) Material Property Table (14, E10.3, 2F10.3, E10.4): One card per 
material type (total NMAT cards) : 



1-4 


Material type number; 


5-14 


Elastic modulus; 


15-24 


Poisson's ratio; 


25-34 


Specific weight; 


35-44 


Coefficient of thermal expansion. 



For plane-strain, reduced values must be used; 

£ - (1 %¥) ^ = (i'-V) <x'zr OC (1 +• y) 

(e) Plotted figure boundary outline array (2014) : For the plotting 

of stress graphs, NUMPB number of corner nodal points which, when con- 
nected outline the structure, or structure subregion, as a series of 
straight lines must be punched, in cyclic order, 20 node numbers per card. 
The starting corner node and direction of travel around the structure, or 



33-34 

35-36 

37-38 

39-40 

41-42 

Notes 
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subregion, boundary are arbitrary. Holes in multiply connected bodies 
cannot be outlined separately. In subregion plots, the elements ( ^ 50) 
outside the plotted figure boundary (skipped elements) must be listed in 
item (1) under Loading. 

(f) Element nodal point array (1014, F10.3): One card per element, 

(total NUMEL cards). 

Cols. 1-4 Element number; 

5-36 Nodal point numbers: 

(I) for quadrilateral: element nodal points in 

counterclockwise order I- J-K-L-M-O-P (Fig. 6.1a). 

The starting corner is arbitrary, except when 
equal stiffnesses are implied (i.e., T1=T in Control 
Card). In this case the starting corners must be 

in the same location for all elements* 

(II) for single triangles: punch nodes I-J-K-L-M-N, 

(Fig. 6.1b), leave cols. 29-36 blank. 

37-40 Element material type, will be set=l if left blank. 

41-50 Element thickness, will be set=] .0 if left blank. 

Note : If a quadrilateral is not convex (not recommended) the entrant 

corner must be either J or K. 

(g) Corner point coordinate array (14, 2F8.3): One card per corner 

point (total NUMCP cards). 

Cols. 1-4 Corner nodal point number; 

5-12 X-coordinate 

13-20 Y -coordinate 

(h) Boundary condition array (214, 2E15«3): One card per restrained 

nodal point (total NUMBC cards). 
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Cols. 1-4 



5-8 



9-23 



24-38 



Nodal point number 

Tag - 0 if point fixed in both directions or initial 
displacement is specified in both the X and Y 
direction then point is fixed. 

= 1 if point is fixed in the X-direction and free 
in the Y-direction or initial displacement is 
specified in the X-direction and point is free 
in the Y-direction. 

= 2 if point is free to move along a line forming 
angle <£) with the X-axis and fixed in a direc- 
tion normal to that line; or if initial dis- 
placement is specified in the Y-direction and 
the point is free in the X-direction. 

Angle in degrees, positive counterclockwise from X-axis 
for type TAG=2 boundary condition. 

Initial X-displacement value for TAG=0 boundary condi- 
tion (for both cases will be set=0.0 if left blank) 
Initial displacement boundary condition value. Y- 
displacement for TAG=0 or 2. X-displacement for TAG=1. 
(will be set = 0.0 if left blank) 



6.2 Loading data 

Each load case must be specified by a data deck initiated by a LOADING 
card; this package follows the structure data deck. A load deck consists 
of the following cards; 

(i) Loading Card (A8) : With the word LOADING punched in Cols. 1-7. 

(j) Title Card (20A4) : Alphameric information in cols. 1-80 to 

identify the load case. 
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Figure 6.1 Element Nodal Point Identification 



Positive Sense Indicated 
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side 2 
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Quadrilateral 
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Figure 6.2 Convention for Element Side Loading 
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(k) Control card (3I4,L2) 
Variable 



Columns 


Name Meaning 


1-4 


NPLD Number of nodal points with concentrated 

forces ; 


5-8 


NELD Number of element sides loaded with dis- 

tributed forces; 


9-12 


NTLD Number of elements undergoing thermal 

increments ; 


13-14 


DENS Logical flag for gravity loading: If a T 



(implies .TRUE,) is punched, gravity forces 
acting along the (-Y) direction are con- 
sidered. 

(1) Stress Contour Graph Indicator Card (714): A positive integer 

punched in any of the first six fields will cause a stress contour plot 
to be produced, on the printer, for the indicated stress component. 



Columns 


Graph 


1-4 


Sigma x 


5-8 


Sigma y 


9-12 


Tau xy 


13-16 


Sigma max 


17-20 


Sigma min 


21-24 


Max shear 



The last field (cols. 25-28) indicates the number NSK 50 of ele- 
ments to be skipped from the plots. If NSK > 0, additional cards must 
follow, specifying the skipped element numbers (2014). Element skipping 
may be used for two different purposes: 

(1) to eliminate small regions of high stress gradients which 
cannot be accurately described by a printer plot; 
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(2) to plot a portion of the structure, which is then amplified. v 
In this case the plotted figure outline array (item e) must 
specify the outline of the subregion. 

(m) Nodal Point Forces (14, 2F8.3): One card per nodal point loaded 

with a concentrated force (no cards if NPLD 35 0). 

Cols. 1-4 Nodal point number; 

5-12 X-load; 

13-20 Y-load. 

(n) Element Side Loads (214, 6F8.3): One card per element side under 

surface traction (no cards if NELD = 0). The convention for positive 
traction and shear is indicated in Figure 6.2. The side variation is 
assumed to be parabolic and specified by the values at points a, b and c 
(in counterclockwise sense). For instance, for side 2 of a quadrilateral: 



corner point 


J, b = midside point N, c = c 


1-4 


Element number; 


5-8 


Side number (see Fig. 6.2); 


9-16 


Normal traction at a (T ) ; 

a 


17-24 


Normal traction at b (T,) ; 

b 


25-32 


Normal traction at c (T ) ; 

c 


33-40 


Surface shear at a; 


41-48 


Surface shear at b; 


49-56 


Surface shear at c. 



These values must be specified per unit of length of the figure and 
per unit of thickness, (eg., PSI) . 

(o) Thermal Increments (14, 4F10.3): One card per element under- 

going temperature changes (no cards if NTLD = 0): 
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Cols. 



1-4 



Element number; 



5-14 Temperature variation at corner I; 

15-24 Temperature variation at corner J; 

25-34 Temperature variation at corner K; 

35-44 Temperature variation at corner L (leave blank for 

a triangle). 

Note: the thermal increment at the centroid of quadrilaterals is assumed 

to be the mean of the corner values, and a linear variation assumed over 
each subtriangle. 

6.3 New problem 

The input of a new problem must follow the last load deck for the 
previous one. For safety, any number of blank cards may be inserted 
before the START card. 
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Recommendations 



The most hearty recommendation is to use the program. It represents 
a powerful engineering tool. 

Recommendations for augmentation of the method are: 

(1) Inclusion of an optional bilinear element stiffness subroutine 
to allow analysis of the ever increasing group of materials that 
have such characteristics. 

(2) Modification to accept linearly varying element thickness. 

(3) A general orthotropic element stiffness subroutine. 
Recommendations for augmentation of the program are: 

(1) A mesh generation package configured for the NPS Computer 
Center IBM Optical Display Unit. 

(2) A contour graph plotting package for use on commercial X-Y 
plotters, and compatible with the program punched card output of 
displacements and stresses. 

(3) Conversion of tape-disk external storage statements in the 
program to direct access statements. 
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APPENDIX 1 



Sample Problem - '’Flounder' 1 plate in tension 



This appendix presents a sample problem analyzed with the program* 

The structure and loading situation are somewhat hypothetical to allow 
presentation of a variety of input data and computed results. Actual 
input data decks and computer output are illustrated. 

A "Flounder" plate in the Naval service is any roughly triangular 
plate, that in its lifetime may inadvertently find itself on the ocean 
floor. The example problem presented in this appendix is the analysis 
of one type of flounder plate. 

The function of the problem plate is that of an attachment member. 

The plate is welded to a "Trolley-Block" apparatus which rides, on cables, 
between ships at sea conducting underway transfer of supplies. The 
flounder plate acts as the attachment point for a small fixture, which in 
turn holds nets or boxes containing the transferred material. The problem 
is the preliminary coarse mesh analysis of the plate. 

Sketches of the complete trolley-block assembly and the loaded 
flounder plate are presented in Figure Al.l. The figure also illustrates 
the idealized form of the plate; nodal point numbers; element numbers; 
nodal coordinates; boundary conditions, and input loading. 

The program input data deck is illustrated in Figure A1.12. Data 
values that are overlined in the figure represent input for which optional 
default values are available. The size of the data fields for input data 
are also illustrated. 

The computed output data is presented in Figure A1.3, where signi- 
ficant items are annotated on the figure. 
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Initial Displ. at Node-5 




FIGURE A l.l SAMPLE PROBLEM STRUCTURE 
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FIGURE AI.2 SAMPLE PROBLEM COMPUTER INPUT DECK 
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FIGURE A 1.3 SAMPLE PROBLEM COMPUTER OUTPUT 
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BOUNCARY CONDITIONS 
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MALONE B0X-3C EXAMPLE PROBLEM - FLOUNDER PLATE UNITS: INCH/KIP/KS 
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NOD/iL FCRCE VECTOR 
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Note : Initial Displacement values appear in the Force 
at their Node and direction of application. See 
column, at Node-5. 



NCCAL PCINT DISPLACEMENTS 
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APPENDIX 2 



Practical Aspects of the Computer Alogrithro 



This appendix presents practical aspects of the computer alogrithm 
for the direct stiffness procedures, equation solver, and the applica- 
tion of boundary conditions. 

The direct stiffness procedure is the process whereby the complete 
structure stiffness matrix is formed from the individual element stiffness 
matrices. The process is conceptually simple, since it requires only the 
systematic addition of the stiffnesses of all the elements in the system. 
Difficulties arise when the concept is transposed to a computer process 
where storage space is at a premium. A firm understanding of the methods 
of assembly and storage of the complete stiffness matrix within the pro- 
gram is of importance for those desiring to modify the existing program. 
These methods control the programming techniques for application of bound- 
ary conditions, and the solution of the equilibrium equations. 

The superposition of each element stiffness, in the formation of the 
complete stiffness matrix, is accomplished by adding its individual terms 
into the complete stiffness matrix according to the nodal point number 
of the elements. This is illustrated in Figure A2.1. 

The program could theoretically construct the complete stiffness 
matrix in the manner illustrated in Fig, A2.1; a two-dimensional array 
could be defined and simple indexing methods utilized to superpose the 
element stiffnesses. The storage requirements for this straightforward 
approach would be a square array with side dimension of (2 x Maximum 
Number of Nodes). The resulting array (2100 x 2100) would occupy the 
total main core storage of 70 computers of the size used for this pro- 
gram, The techniques used to alleviate the storage problem utilizes 
three characteristics of the complete stiffness matrix. 

The complete stiffness matrix is 
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FIG. A2.I DIRECT STIFFNESS PROCEDURE 
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(1) Symmetric 

(2) Banded 

(3) Sparsely populated 

These properties are illustrated in Figure A2.2. 




Figure A2.2 Complete Stiffness Matrix Characteristics 
Symmetry permits a reduction of approximately one-half in the storage 
required, since only a triangular half of the matrix is required to 
retain significant data. The banded and sparseness properties allow 
further savings, since only data between the diagonal and the band limits 
is non-zero and required in computation. The cross-hatched portion of 
Fig. A2.2 shows the actual quantity of storage required. The storage space 
problem is only partially solved, since space approximately equal to an 
array size dimensioned (% Bandwidth x Diagonal Length) is still required*. 

In the program the maximum ha If -bandwidth is 160 and the maximum diagonal 
length (equals a side dimension) is 2100. The resulting (160 x 2100) ar- 
ray is small enough to be stored on an external disk storage unit, but 
cannot be housed completely in main core (5 would be required) where high 
speed calculations, independent of internal/external storage transfer 
time, can be performed. 
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Considerations of the equation solving method and transfer time 
reduction dictate the final techniques used in the storage and manipula- 
tion of the data contained in the complete stiffness matrix. These two 
considerations will be discussed before proceeding. 

The equation solver used in the program employs Cholesky's algorithm 
[4], [10]. The fundamental concept of the solution process is to perform 
operations on the complete stiffness matrix to reduce it to triangular 
form, after which the unknown displacements are found by back substitu- 
tion procedures. The limited coupling characteristic of the set of 
equations allows the reduction in form to be accomplished in a manner 
where a subgrouping of equations can be reduced independently in a pro- 
cess termed triangular decomposition. The decomposition process does 
not modify the load vector in the equilibrium equations, and once com- 
pleted, the decomposed form of the complete stiffness matrix may be 
stored for possible use with additional load cases. This accounts for 
the greatly reduced running time for additional load cases. The sub- 
groups of equations are handled as individual blocks of equations. 

Figure A2.3 shows a block structure used for decomposing the complete 
stiffness matrix. The figure also illustrates two points concerning 
block structure; (1) element stiffness contributions may enter two dif- 
ferent blocks, and (2) in order to maintain equal block size, the first 
must contain a zero value area. 

The transfer time between external and internal core storage is a 
function of the data address determination time and the physical time 
required to transfer the data. Where many individual units or small 
groups of data are being moved from within a larger data set (the situa- 
tion in our case), the address computation time is usually the largest 
contributor to the overall time. To reduce address time to a minimum, 
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zero value area 



Block 

Numbers 
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A single element may 




Figure A2.3 Block Structure used for decomposing [K] 
the use of singly dimensioned arrays, with their attendent short address 
calculations, is required. 

Returning to the complete stiffness matrix, with its maximum data 
storage requirement of a (160 x 2100) array, the considerations of the 
preceeding paragraphs are applied. The final storage technique for the 
resulting matrix data calls for data in block form with each block stored 
individually on an external disk unit. Each block is transferred, as 
required, to main core storage for the formation, modification, or solu- 
tion procedures. Each block is stored in a one-dimensional array form. 
The overall storage method is illustrated in Figure A2.4, where an orig- 
inal complete stiffness matrix is shown, along with a virtual array 
depicting the actual data that requires storage, and the final one- 
dimensional array form of one of the stored blocks. Block divisions are 
indicated in the figure. The rows of the virtual array, and the order of 
data in the one -dimensional array, are the columns (read bottom to top) 
of the upper half -bandwidth of the complete stiffness matrix. The first 
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FIG. A2.4 COMPLETE STIFFNESS MATRIX STORAGE METHOD 
AND BOUNDARY CONDITION APPLICATION 



column of the virtual array is the main diagonal of the complete matrix. 
The original rows become diagonals in the virtual array. 

The application of boundary conditions to the structure equilibrium 
equations involves adjustment of the complete stiffness matrix to re- 
flect the constrained nodes or initially displaced and constrained nodes. 
In the case of constrained nodes, the condition is accounted for by strik- 
ing out the rows and columns corresponding to the degrees of freedom 
associated with the boundary condition and replacing the corresponding 
diagonal terms with a non-zero value. If initial displacements are speci- 
fied the corresponding column of the matrix must be multiplied by the 
initial displacement value, and the resulting vector subtracted from the 
load vector before the above constraint procedures are carried out. In 
addition, the known displacement value must be inserted into the displace- 
ment vector. The boundary condition application procedure is illustrated 
in Figure A2.4, where a 5-element bar, fixed at the ends (nodes 3 and 28); 
with an initial positive Y-displacement at node 13, is loaded with a 
negative Y-direction concentrated force, P, at node 16. 

The carry over of boundary conditions into the virtual storage array 
and one block of the final storage form, is also illustrated in Figure 
A2.4. Since the stiffness data is actually stored in a singly dimensioned 
array, special indexing and address accounting procedures (see subroutine) 
FORMK) must be used to allow complete row and column boundary condition 
procedures to be applied in situations where two data blocks are involved; 
each block being available for manipulation only when in main core. 
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Appendix 3 



Program (PSELST) Listing 
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SUBROUTINE SETUP 
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THIS SUBROUTINE ASSEMBLES AND CONDENSES THE STIFFNESS MATRIX OF 
A CUADR ILATERAL FORMED BY FOUR LST ELEMENTS 

************************** 4 ********* ****** ********************* 
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CONDENSATION OF INTERNAL NOOAL POINTS 
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